#include <iterator>
#include <vector>
#include <algorithm>

/* 
   Compares the *first values of two pairs using comp.
 */
template < typename ForwardIterator, typename StrictWeakOrdering >
class __multiway_merge_compare {
    StrictWeakOrdering comp;
  public:
    __multiway_merge_compare( StrictWeakOrdering &_comp ) : comp(_comp) {}
    bool operator()( std::pair< ForwardIterator, ForwardIterator > const &a,
		     std::pair< ForwardIterator, ForwardIterator > const &b ) {
	// push empty iterator pairs to end.
	if ( a.first == a.second ) { 
	    return b.first != b.second;
	}
	if ( b.first == b.second ) {
	    return false;
	}
	return comp( *(a.first), *(b.first) );
    }
};

template < typename InputIteratorIterator,
	   typename OutputIterator,
	   typename StrictWeakOrdering >
void multiway_merge( InputIteratorIterator begins, 
		     InputIteratorIterator begins_end,
		     InputIteratorIterator ends,
		     OutputIterator out,
		     StrictWeakOrdering comp ) 
{
    if ( begins == begins_end ) return;
    typedef typename std::iterator_traits<InputIteratorIterator>::value_type type_iterator;
    typedef std::pair< type_iterator, type_iterator >                        pair;
    typedef std::vector< pair >                                              pair_vector;
    pair_vector heap;
    std::transform( begins, begins_end, ends,
    		    std::back_inserter( heap ),
		    std::make_pair<type_iterator,type_iterator> );
    
    __multiway_merge_compare<type_iterator, StrictWeakOrdering> lcomp(comp);
    std::make_heap( heap.begin(), heap.end(), lcomp );

    typedef typename pair_vector::iterator::difference_type diff_t;
    typedef typename pair_vector::iterator::value_type value_t;
    const diff_t last( heap.end() - heap.begin() );
    while( heap.front().first != heap.front().second ) {
	*out++ = *(heap.front().first++);
	// sift down
	const value_t tmp(heap.front());
	diff_t cur = 0;
	while ( 2*cur+1 < last ) {
	    diff_t next = 2*cur+1;
	    if ( next+1 < last && lcomp(heap[next], heap[next+1])) {
		next++;
	    }
	    if ( lcomp( heap[next], tmp ) ) {
		break;
	    }
	    heap[cur] = heap[next];
	    cur = next;
	}
	if ( cur ) heap[cur] = tmp;
    }
}
